options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(1,4))
      drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
      drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
      drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
      drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}

stan basic

normal distribution

ex1.stan

normal distribution

data{
  int N;
  vector[N] y;
}

parameters{
  real m;
  real<lower=0> s;
}

model{
  y~normal(m,s);
}
mdl=cmdstan_model('./ex1.stan')
y=rnorm(10,2,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   -0.42    1.79    2.29    2.17    3.16    3.86
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')  

data=list(N=length(y),y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -178.954 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       15      -7.14868    0.00640488   0.000400752           1           1       17    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.3 seconds.
mle
##  variable estimate
##      lp__    -7.15
##      m        2.17
##      s        1.24
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
mcmc
##  variable  mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##      lp__ -8.05  -7.71 1.12 0.83 -10.35 -6.97 1.00      746      687
##      m     2.17   2.18 0.52 0.47   1.30  2.99 1.00      647      609
##      s     1.55   1.47 0.46 0.38   0.99  2.35 1.00     1154     1077
mcmc$metadata()$stan_variables
## [1] "lp__" "m"    "s"
mcmc$metadata()$model_params
## [1] "lp__" "m"    "s"
seeMCMC(mcmc)
##  variable  mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##      lp__ -8.05  -7.71 1.12 0.83 -10.35 -6.97 1.00      746      687
##      m     2.17   2.18 0.52 0.47   1.30  2.99 1.00      647      609
##      s     1.55   1.47 0.46 0.38   0.99  2.35 1.00     1154     1077

poisson distribution

y=rpois(20,3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    2.00    4.00    3.40    4.25    8.00
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=length(y),y=y)

ex2-1.stan

in case fit poisson dist. to normal dist.

data{
  int N;
  vector[N] y;
}

parameters{
  real<lower=0> m;
  real<lower=0> s;
}

model{
  y~normal(m,s);
}

generated quantities{
  vector[N] y1;
  for(i in 1:N)
    y1[i]=normal_rng(m,s);
}
# when fitting poisson dist. sample to normal dist.
mdl=cmdstan_model('./ex2-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -29.1211 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5      -22.0597    0.00121336   0.000370207           1           1       10    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -22.06
##    m          3.40
##    s          1.83
##    y1[1]      3.25
##    y1[2]      0.83
##    y1[3]      4.25
##    y1[4]      5.06
##    y1[5]      1.06
##    y1[6]      3.24
##    y1[7]      5.54
##    y1[8]      1.40
##    y1[9]      5.38
##    y1[10]     6.61
##    y1[11]     1.67
##    y1[12]     4.02
##    y1[13]     2.19
##    y1[14]     7.61
##    y1[15]     5.49
##    y1[16]     4.65
##    y1[17]     6.51
##    y1[18]     4.31
##    y1[19]     5.76
##    y1[20]     1.75
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -21.28 -20.98 1.04 0.75 -23.47 -20.27 1.01      704     1032
##    m        3.41   3.40 0.45 0.43   2.70   4.13 1.00     1628     1102
##    s        2.02   1.98 0.37 0.35   1.52   2.68 1.00      796     1104
##    y1[1]    3.29   3.29 2.16 2.05  -0.25   6.80 1.00     1889     1794
##    y1[2]    3.32   3.32 2.14 2.11  -0.22   6.80 1.00     1879     1817
##    y1[3]    3.41   3.37 2.10 2.00   0.03   6.89 1.00     1983     1931
##    y1[4]    3.37   3.39 2.12 2.03  -0.18   6.93 1.00     2204     2012
##    y1[5]    3.43   3.47 2.06 1.98   0.08   6.73 1.00     1699     1787
##    y1[6]    3.48   3.53 2.10 1.99  -0.01   6.82 1.00     1951     1956
##    y1[7]    3.36   3.42 2.10 1.97  -0.22   6.79 1.00     2156     1949
##    y1[8]    3.37   3.35 2.12 2.00  -0.11   6.90 1.00     1963     1816
##    y1[9]    3.35   3.37 2.06 1.96  -0.05   6.70 1.00     1978     1646
##    y1[10]   3.41   3.34 2.13 2.01   0.06   6.85 1.00     1925     1975
##    y1[11]   3.51   3.48 2.14 2.02   0.06   7.06 1.00     1798     1532
##    y1[12]   3.38   3.37 2.11 2.06  -0.11   6.90 1.00     2274     2035
##    y1[13]   3.44   3.47 2.11 2.03  -0.09   6.98 1.00     2130     1933
##    y1[14]   3.38   3.37 2.07 2.07  -0.06   6.58 1.00     1823     1712
##    y1[15]   3.47   3.48 2.15 2.09  -0.19   6.96 1.00     2006     1840
##    y1[16]   3.42   3.46 2.12 1.99  -0.04   6.82 1.00     1887     1858
##    y1[17]   3.42   3.40 2.10 1.95  -0.04   6.81 1.00     1896     1894
##    y1[18]   3.41   3.44 2.12 2.04  -0.13   6.79 1.00     1867     1741
##    y1[19]   3.44   3.39 2.05 2.01   0.21   6.85 1.00     2161     1916
##    y1[20]   3.38   3.36 2.18 2.13  -0.26   7.06 1.00     1930     2011

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1')     

ex2-2.stan

fit to poisson dist.

data{
  int N;
  array[N] int y;
}

parameters{
  real<lower=0> l;
}

model{
  y~poisson(l);
}

generated quantities{
  array[N] int y1;
  for(i in 1:N)
    y1[i]=poisson_rng(l);
}
mdl=cmdstan_model('./ex2-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -4.6355 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6       15.2167   8.35723e-05   1.05621e-05           1           1        7    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__      15.22
##    l          3.40
##    y1[1]      4.00
##    y1[2]      6.00
##    y1[3]      2.00
##    y1[4]      4.00
##    y1[5]      3.00
##    y1[6]      8.00
##    y1[7]      3.00
##    y1[8]      4.00
##    y1[9]      4.00
##    y1[10]     3.00
##    y1[11]     3.00
##    y1[12]     3.00
##    y1[13]     6.00
##    y1[14]     0.00
##    y1[15]     5.00
##    y1[16]     2.00
##    y1[17]     4.00
##    y1[18]     2.00
##    y1[19]     6.00
##    y1[20]     2.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##    lp__   15.96  16.24 0.72 0.29 14.56 16.45 1.00      852     1068
##    l       3.43   3.42 0.41 0.39  2.79  4.15 1.00      731      939
##    y1[1]   3.52   3.00 1.92 1.48  1.00  7.00 1.00     1689     1711
##    y1[2]   3.44   3.00 1.86 1.48  1.00  7.00 1.00     1881     1926
##    y1[3]   3.45   3.00 1.86 1.48  1.00  7.00 1.00     1740     1832
##    y1[4]   3.40   3.00 1.91 1.48  1.00  7.00 1.00     1919     1988
##    y1[5]   3.48   3.00 1.88 1.48  1.00  7.00 1.00     1692     1853
##    y1[6]   3.44   3.00 1.91 1.48  1.00  7.00 1.00     2017     1814
##    y1[7]   3.50   3.00 1.91 1.48  1.00  7.00 1.00     1716     2033
##    y1[8]   3.38   3.00 1.89 1.48  1.00  7.00 1.00     1744     1959
##    y1[9]   3.47   3.00 1.93 1.48  1.00  7.00 1.00     1927     2035
##    y1[10]  3.38   3.00 1.91 1.48  1.00  7.00 1.00     1710     1892
##    y1[11]  3.39   3.00 1.87 1.48  1.00  7.00 1.00     1799     1831
##    y1[12]  3.45   3.00 1.91 1.48  1.00  7.00 1.00     1892     1863
##    y1[13]  3.42   3.00 1.89 1.48  1.00  7.00 1.00     1838     1913
##    y1[14]  3.44   3.00 1.90 1.48  1.00  7.00 1.00     1816     1846
##    y1[15]  3.41   3.00 1.92 1.48  1.00  7.00 1.00     1684     1636
##    y1[16]  3.33   3.00 1.89 1.48  1.00  7.00 1.00     1805     1880
##    y1[17]  3.42   3.00 1.88 1.48  1.00  7.00 1.00     1893     1826
##    y1[18]  3.35   3.00 1.84 1.48  1.00  7.00 1.00     1740     1905
##    y1[19]  3.46   3.00 1.85 1.48  1.00  7.00 1.00     1868     1917
##    y1[20]  3.40   3.00 1.90 1.48  1.00  7.00 1.00     1633     2023

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

difference test

ex3.stan

mean difference of 2 normal distributions

data{
  int N;
  vector[N] a;
  vector[N] b;
}
parameters{
  real ma;
  real<lower=0> sa;
  real mb;
  real<lower=0> sb;
}
model{
  a~normal(ma,sa);
  b~normal(mb,sb);
}
generated quantities{
  real d;
  d=mb-ma;
}
n=20
a=rnorm(n,10,1)
b=rnorm(n,11,2)
data=list(N=n,a=a,b=b)

mdl=cmdstan_model('./ex3.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -1244.97 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       26      -24.8565   0.000370989   0.000220677           1           1       43    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -24.86
##      ma       9.89
##      sa       0.69
##      mb      10.57
##      sb       1.84
##      d        0.68
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -26.71 -26.36 1.51 1.35 -29.77 -24.93 1.01      796     1389
##      ma     9.89   9.89 0.18 0.17   9.60  10.19 1.00     2099     1257
##      sa     0.76   0.74 0.14 0.13   0.57   1.01 1.00     2129     1503
##      mb    10.56  10.56 0.45 0.46   9.81  11.29 1.01      644      596
##      sb     2.03   1.99 0.35 0.33   1.54   2.71 1.00     2386     1454
##      d      0.66   0.65 0.49 0.49  -0.14   1.50 1.01      727      746

d=mcmc$draws('d')
d
## # A draws_array: 500 iterations, 4 chains, and 1 variables
## , , variable = d
## 
##          chain
## iteration     1     2    3    4
##         1  0.42 -0.15 0.81 0.67
##         2  0.87  1.48 1.15 0.54
##         3 -0.11  1.04 1.30 1.22
##         4  0.89  0.66 1.50 0.58
##         5  0.79  0.98 1.34 0.49
## 
## # ... with 495 more iterations
mean(d>0)
## [1] 0.911

normal regression

ex4-1.stan

single normal regression

data{
  int N;
  vector[N] x;
  vector[N] y;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}
# make prediction with R
n=20
x=runif(n,0,100)
y=rnorm(n,x*3+10,10)
par(mfrow=c(1,1))
plot(x,y)

data=list(N=n,x=x,y=y)

mdl=cmdstan_model('./ex4-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -1.72607e+06 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       42       -52.917    0.00706613   0.000642451      0.9445     0.09445      117    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -52.92
##      b0      14.59
##      b1       3.00
##      s        8.55
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -52.44 -52.09 1.37 1.16 -55.13 -50.96 1.01      320      508
##      b0    14.46  14.62 4.57 3.94   7.08  22.11 1.01      219      167
##      b1     3.00   3.00 0.08 0.07   2.87   3.13 1.01      301      398
##      s      9.75   9.53 1.74 1.64   7.40  13.02 1.00      901      892

b0=mcmc$draws('b0')
b1=mcmc$draws('b1')
b=tibble(b0=c(b0),b1=c(b1)) |> sample_n(100)
x1=runif(nrow(b),min(x),max(x))
xy=tibble(x=x1,y=b$b0+b$b1*x1)
par(mfrow=c(1,1))
plot(xy)

ex4-2.stan

single normal regression, prediction from new explanatory variable

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}

generated quantities{
  vector[N1] m;
  vector[N1] y1;
  for(i in 1:N1){
    m[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m[i],s);
  }
}
# make prediction with stan
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex4-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -1.88247e+06 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       48       -52.917   0.000245332    0.00357952      0.4582      0.4582      107    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -52.92
##    b0        14.59
##    b1         3.00
##    s          8.55
##    m[1]      28.46
##    m[2]      58.37
##    m[3]      88.28
##    m[4]     118.19
##    m[5]     148.10
##    m[6]     178.01
##    m[7]     207.93
##    m[8]     237.84
##    m[9]     267.75
##    m[10]    297.66
##    y1[1]     37.39
##    y1[2]     61.60
##    y1[3]     85.19
##    y1[4]    125.25
##    y1[5]    150.33
##    y1[6]    168.23
##    y1[7]    204.06
##    y1[8]    236.88
##    y1[9]    285.83
##    y1[10]   297.70
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 0.3 seconds.
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 0.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "s"    "m"    "y1"
seeMCMC(mcmc,'m',ch=F)
##  variable   mean median    sd   mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -52.32 -52.00  1.27  0.99 -54.78 -50.94 1.00      580      942
##    b0      14.89  14.65  4.15  4.17   8.57  21.85 1.01      226      267
##    b1       2.99   3.00  0.07  0.07   2.87   3.11 1.01      303      494
##    s        9.71   9.45  1.72  1.58   7.38  12.86 1.00     1038      955
##    m[1]    28.75  28.48  3.86  3.86  22.80  35.18 1.01      226      278
##    m[2]    58.63  58.43  3.27  3.31  53.45  64.03 1.01      232      303
##    m[3]    88.52  88.42  2.75  2.69  84.10  93.07 1.01      253      367
##    m[4]   118.40 118.35  2.36  2.30 114.58 122.34 1.01      317      536
##    m[5]   148.29 148.29  2.16  2.07 144.70 151.76 1.00      532      879
##    m[6]   178.17 178.21  2.21  2.10 174.57 181.63 1.00     1326     1338
##    m[7]   208.06 208.04  2.50  2.37 204.00 212.02 1.00     2336     1609
##    m[8]   237.94 237.94  2.95  2.75 233.29 242.69 1.00     1851     1428
##    m[9]   267.83 267.79  3.50  3.22 262.33 273.36 1.00     1162     1316
##    m[10]  297.72 297.70  4.11  3.75 291.17 304.27 1.00      855     1183
##    y1[1]   28.45  28.26 10.42  9.93  11.31  46.14 1.00     1016     1737
##    y1[2]   58.46  58.79 10.41 10.21  41.27  75.00 1.00     1099     1712
##    y1[3]   88.45  88.56 10.27  9.93  71.27 104.58 1.00     1455     1766
##    y1[4]  118.69 118.58  9.93  9.56 102.81 135.43 1.00     1561     1924
##    y1[5]  148.46 148.51  9.81  9.52 132.63 164.67 1.00     1816     1728
##    y1[6]  178.02 178.02 10.07  9.49 161.34 194.70 1.00     1993     2016
##    y1[7]  207.81 207.74 10.06  9.85 191.13 224.36 1.00     2058     1725
##    y1[8]  237.68 237.82 10.40 10.19 220.94 255.18 1.00     2166     2010
##    y1[9]  268.08 268.11 10.63 10.46 250.08 285.64 1.00     1634     2101
##    y1[10] 297.88 298.02 10.70 10.19 280.66 315.03 1.00     1519     1784

m=mcmc$draws('m')
summary(m)
## # A tibble: 10 × 10
##    variable  mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
##    <chr>    <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
##  1 m[1]      28.7   28.5  3.86  3.86  22.8  35.2 1.01      227.     279.
##  2 m[2]      58.6   58.4  3.27  3.31  53.5  64.0 1.01      233.     304.
##  3 m[3]      88.5   88.4  2.75  2.69  84.1  93.1 1.01      253.     367.
##  4 m[4]     118.   118.   2.36  2.30 115.  122.  1.01      318.     536.
##  5 m[5]     148.   148.   2.16  2.07 145.  152.  1.00      533.     879.
##  6 m[6]     178.   178.   2.21  2.10 175.  182.  0.999    1327.    1338.
##  7 m[7]     208.   208.   2.50  2.37 204.  212.  1.00     2336.    1609.
##  8 m[8]     238.   238.   2.95  2.75 233.  243.  1.00     1851.    1429.
##  9 m[9]     268.   268.   3.50  3.22 262.  273.  1.00     1162.    1317.
## 10 m[10]    298.   298.   4.11  3.75 291.  304.  1.00      855.    1184.
smm=summary(m)

y1=mcmc$draws('y1')
summary(y1)
## # A tibble: 10 × 10
##    variable  mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
##    <chr>    <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
##  1 y1[1]     28.5   28.3 10.4   9.93  11.3  46.1 1.00     1017.    1738.
##  2 y1[2]     58.5   58.8 10.4  10.2   41.3  75.0 1.00     1099.    1712.
##  3 y1[3]     88.5   88.6 10.3   9.93  71.3 105.  1.00     1456.    1766.
##  4 y1[4]    119.   119.   9.93  9.56 103.  135.  1.00     1562.    1924.
##  5 y1[5]    148.   149.   9.81  9.52 133.  165.  1.00     1816.    1728.
##  6 y1[6]    178.   178.  10.1   9.49 161.  195.  1.00     1994.    2017.
##  7 y1[7]    208.   208.  10.1   9.85 191.  224.  1.00     2059.    1725.
##  8 y1[8]    238.   238.  10.4  10.2  221.  255.  1.00     2166.    2011.
##  9 y1[9]    268.   268.  10.6  10.5  250.  286.  0.999    1635.    2102.
## 10 y1[10]   298.   298.  10.7  10.2  281.  315.  1.00     1519.    1785.
smy=summary(y1)

xy=tibble(x=x1,m=smm$median,ml5=smm$q5,mu5=smm$q95,yl5=smy$q5,yu5=smy$q95)

par(mfrow=c(1,1))
xlim=c(min(x),max(x))
ylim=c(min(y),max(y))
plot(x,y,
     xlim=xlim,ylim=ylim, xlab='x',ylab='y',col='red')
par(new=T)
plot(xy$x,xy$m,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='black')
par(new=T)
plot(xy$x,xy$ml5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$mu5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$yl5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')
par(new=T)
plot(xy$x,xy$yu5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=m),xy,col='black')+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')

ex4-3.stan

single normal regression, prediction from original explanatory variable

data{
  int N;
  vector[N] x;
  vector[N] y;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}

generated quantities{
  vector[N] m;
  vector[N] y1;
  for(i in 1:N){
    m[i]=b0+b1*x[i];
    y1[i]=normal_rng(m[i],s);
  }
}
data=list(N=n,x=x,y=y)

mdl=cmdstan_model('./ex4-3.stan')

mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 0.2 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,'m',ch=F)
##  variable   mean median    sd   mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -52.25 -51.94  1.20  0.94 -54.64 -50.95 1.00      509      773
##    b0      14.77  14.68  4.03  4.09   8.86  21.52 1.01      177      193
##    b1       2.99   2.99  0.07  0.07   2.87   3.10 1.01      244      381
##    s        9.61   9.39  1.70  1.60   7.30  12.79 1.00     1039     1025
##    m[1]    51.16  51.06  3.35  3.42  46.28  56.72 1.01      182      192
##    m[2]   267.85 267.91  3.21  3.12 262.49 273.13 1.00     1377     1216
##    m[3]   137.28 137.24  2.21  2.17 133.63 140.89 1.00      337      676
##    m[4]   227.00 227.02  2.58  2.52 222.79 231.27 1.00     2040     1404
##    m[5]   275.20 275.27  3.34  3.21 269.59 280.67 1.00     1251     1229
##    m[6]    74.78  74.65  2.95  3.06  70.39  79.58 1.01      192      205
##    m[7]   261.01 261.07  3.10  2.99 255.83 266.09 1.00     1518     1294
##    m[8]    28.63  28.51  3.77  3.86  23.19  34.92 1.01      178      194
##    m[9]   169.78 169.78  2.14  2.05 166.37 173.32 1.00      773     1077
##    m[10]  189.54 189.57  2.21  2.19 185.95 193.19 1.00     1493     1417
##    m[11]   48.19  48.08  3.41  3.47  43.24  53.84 1.01      181      199
##    m[12]  290.13 290.19  3.61  3.43 284.05 296.04 1.00     1051     1117
##    m[13]  181.12 181.13  2.17  2.11 177.66 184.72 1.00     1098     1372
##    m[14]  297.48 297.56  3.75  3.57 291.21 303.56 1.00      975     1067
##    m[15]  143.22 143.19  2.18  2.13 139.67 146.80 1.00      383      708
##    m[16]  160.10 160.09  2.13  2.07 156.62 163.61 1.00      592     1039
##    m[17]   65.63  65.51  3.10  3.21  61.05  70.70 1.01      187      190
##    m[18]  224.37 224.37  2.55  2.51 220.26 228.61 1.00     2057     1363
##    m[19]   77.34  77.19  2.91  3.03  72.98  82.02 1.01      194      220
##    m[20]   71.42  71.29  3.01  3.10  66.91  76.28 1.01      190      202
##    y1[1]   50.94  50.71 10.43 10.27  34.30  68.01 1.00     1178     1556
##    y1[2]  267.63 267.83 10.16  9.59 250.63 284.15 1.00     2146     1856
##    y1[3]  137.27 137.50 10.22  9.91 120.33 153.66 1.00     1584     1858
##    y1[4]  227.16 227.24 10.37 10.20 209.93 244.14 1.00     2000     1772
##    y1[5]  274.90 274.96 10.30 10.17 258.01 291.44 1.00     1715     1649
##    y1[6]   75.05  74.89 10.17 10.18  58.26  91.83 1.00     1471     1894
##    y1[7]  261.09 261.15  9.94  9.35 244.81 277.62 1.00     1921     1956
##    y1[8]   28.65  28.41 10.65 10.02  11.43  46.09 1.00     1028     1740
##    y1[9]  169.62 169.62 10.08  9.98 153.09 185.73 1.00     1849     1731
##    y1[10] 189.47 189.50  9.96  9.42 172.91 205.61 1.00     1742     1931
##    y1[11]  48.17  48.31 10.19  9.95  30.87  64.46 1.00     1214     1608
##    y1[12] 289.68 289.67 10.03  9.51 273.53 306.54 1.00     1938     1805
##    y1[13] 181.18 181.20 10.08  9.32 164.58 197.39 1.00     2188     2001
##    y1[14] 297.13 297.36 10.34  9.90 280.00 314.44 1.00     1688     1711
##    y1[15] 143.04 143.08 10.12  9.68 127.03 159.71 1.00     1975     2046
##    y1[16] 160.37 160.53 10.12  9.51 143.45 176.60 1.00     1966     1857
##    y1[17]  65.59  65.74 10.09  9.82  49.02  82.23 1.00     1011     1431
##    y1[18] 224.05 223.88 10.30  9.99 207.18 241.06 1.00     1981     1523
##    y1[19]  77.26  77.59 10.33  9.60  60.38  93.23 1.00     1203     1778
##    y1[20]  71.77  71.56 10.21 10.16  55.40  88.60 1.00     1503     1567

y1=mcmc$draws('y1')
smy=summary(y1)

par(mfrow=c(1,1))
plot(y,smy$median,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(y,smy$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(y-smy$median,xlab='obs.-prd.',main='residual')
density(y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

estimate correlation

ex4-41.stan

use covariance matrix and multinormal distribution

data{
  int N;
  array[N] vector[2] y;
}

parameters{
  vector[2] m;
  vector<lower=0>[2] s;
  real<lower=-1,upper=1> r;
}

transformed parameters{
  cov_matrix[2] cv;
  cv[1,1]=s[1]^2;
  cv[2,2]=s[2]^2;
  cv[1,2]=r*s[1]*s[2];
  cv[2,1]=r*s[1]*s[2];
}

model{
  y~multi_normal(m,cv);
}

ex4-42.stan

use covariance matrix and wishart distribution

data{
  int N;
  matrix[2,2] dp;
}

parameters{
  vector<lower=0>[2] s;
  real<lower=-1,upper=1> r;
}

transformed parameters{
  cov_matrix[2] cv;
  cv[1,1]=s[1]^2;
  cv[2,2]=s[2]^2;
  cv[1,2]=r*s[1]*s[2];
  cv[2,1]=r*s[1]*s[2];
}

model{
  dp~wishart(N-1,cv);
}

ex4-43.stan

use correlation matrix and wishart distribution

data{
  int N;
  int M;
  matrix[M,M] dp;
}

parameters{
  vector<lower=0>[M] s;
  corr_matrix[M] r;
}

transformed parameters{
  cov_matrix[M] cv;
  cv=quad_form_diag(r,s);
}

model{
  dp~wishart(N-1,cv);
}

ex4-44.stan use covariance matrix and multinormal distribution

data{
  int N;
  int K;
  array[N] vector[K] y;
}
parameters{
  vector[K] m;
  cov_matrix[K] cov;
}
model{
  y~multi_normal(m,cov);
}
library(MASS)
n=20
y=mvrnorm(n,c(0,0),matrix(c(1,0.7,0.7,1),2),2)
cov(y)
##       [,1]  [,2]
## [1,] 0.798 0.626
## [2,] 0.626 0.770
cor(y)
##       [,1]  [,2]
## [1,] 1.000 0.798
## [2,] 0.798 1.000
plot(y)

#estimate covariance
mdl=cmdstan_model('./ex4-41.stan')

data=list(N=n,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -77.4925 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       26      -3.96424   1.35683e-05   0.000643409      0.4565      0.4565       30    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__       -3.96
##   m[1]        0.08
##   m[2]       -0.34
##   s[1]        0.87
##   s[2]        0.86
##   r           0.80
##   cv[1,1]     0.76
##   cv[2,1]     0.59
##   cv[1,2]     0.59
##   cv[2,2]     0.73
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable  mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##   lp__    -8.65  -8.31 1.71 1.52 -11.86 -6.54 1.00      661      864
##   m[1]     0.09   0.09 0.23 0.22  -0.29  0.44 1.00     1057      908
##   m[2]    -0.33  -0.33 0.22 0.21  -0.69  0.02 1.00      985     1080
##   s[1]     0.96   0.94 0.16 0.15   0.73  1.25 1.00     1002      926
##   s[2]     0.94   0.92 0.17 0.15   0.71  1.25 1.00      957      951
##   r        0.76   0.77 0.11 0.09   0.55  0.89 1.01      475      577
##   cv[1,1]  0.94   0.88 0.34 0.28   0.53  1.57 1.00     1002      926
##   cv[2,1]  0.71   0.65 0.30 0.24   0.34  1.26 1.00      600      647
##   cv[1,2]  0.71   0.65 0.30 0.24   0.34  1.26 1.00      600      647
##   cv[2,2]  0.91   0.84 0.35 0.27   0.51  1.57 1.00      957      951

#estimate correlation,covariance
#deviation product sum matirix folows Wishart dist.
#deviation product sum = covariance * n-1 or n

mdl=cmdstan_model('./ex4-42.stan')

data=list(N=n,dp=cov(y)*(n-1))

mle=mdl$optimize(data=data)
## Initial log joint probability = -62.2028 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       16       -4.7406   3.36999e-05   4.25106e-05           1           1       21    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__       -4.74
##   s[1]        0.89
##   s[2]        0.88
##   r           0.80
##   cv[1,1]     0.80
##   cv[2,1]     0.63
##   cv[1,2]     0.63
##   cv[2,2]     0.77
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable  mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##   lp__    -8.20  -7.85 1.27 1.05 -10.71 -6.84 1.00      708      983
##   s[1]     0.95   0.93 0.17 0.15   0.72  1.28 1.00      806      862
##   s[2]     0.93   0.91 0.17 0.15   0.71  1.25 1.00      825     1017
##   r        0.75   0.77 0.11 0.10   0.54  0.89 1.01      406      475
##   cv[1,1]  0.93   0.86 0.34 0.27   0.53  1.64 1.00      806      862
##   cv[2,1]  0.69   0.63 0.30 0.24   0.34  1.29 1.00      529      607
##   cv[1,2]  0.69   0.63 0.30 0.24   0.34  1.29 1.00      529      607
##   cv[2,2]  0.90   0.82 0.34 0.27   0.51  1.56 1.00      825     1017

#estimate correlation,covariance
mdl=cmdstan_model('./ex4-43.stan')

m=2
data=list(N=n,M=m,dp=cov(y)*(n-1))

mle=mdl$optimize(data=data)
## Initial log joint probability = -36.6064 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       12       -4.7406   0.000412103   1.97173e-05           1           1       14    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__       -4.74
##   s[1]        0.89
##   s[2]        0.88
##   r[1,1]      1.00
##   r[2,1]      0.80
##   r[1,2]      0.80
##   r[2,2]      1.00
##   cv[1,1]     0.80
##   cv[2,1]     0.63
##   cv[1,2]     0.63
##   cv[2,2]     0.77
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable  mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##   lp__    -7.51  -7.17 1.27 1.05 -10.01 -6.15 1.01      636      858
##   s[1]     0.95   0.93 0.17 0.15   0.72  1.26 1.01     1065      934
##   s[2]     0.93   0.91 0.17 0.15   0.72  1.24 1.00     1040      715
##   r[1,1]   1.00   1.00 0.00 0.00   1.00  1.00   NA       NA       NA
##   r[2,1]   0.75   0.77 0.11 0.10   0.54  0.89 1.01      793      530
##   r[1,2]   0.75   0.77 0.11 0.10   0.54  0.89 1.01      793      530
##   r[2,2]   1.00   1.00 0.00 0.00   1.00  1.00   NA       NA       NA
##   cv[1,1]  0.93   0.87 0.34 0.28   0.52  1.59 1.01     1065      934
##   cv[2,1]  0.69   0.63 0.31 0.24   0.33  1.26 1.00      805      760
##   cv[1,2]  0.69   0.63 0.31 0.24   0.33  1.26 1.00      805      760
##   cv[2,2]  0.90   0.83 0.35 0.26   0.52  1.54 1.00     1040      715

mdl=cmdstan_model('./ex4-44.stan')

k=2
data=list(N=n,K=k,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -261.095 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       20      -3.96424   7.74109e-06   0.000432026      0.3788      0.3788       23    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##  lp__        -3.96
##  m[1]         0.08
##  m[2]        -0.34
##  cov[1,1]     0.76
##  cov[2,1]     0.59
##  cov[1,2]     0.59
##  cov[2,2]     0.73
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable  mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##  lp__     -6.97  -6.60 1.84 1.49 -10.54 -4.79 1.01      706      753
##  m[1]      0.07   0.07 0.24 0.23  -0.32  0.46 1.00      926      815
##  m[2]     -0.34  -0.34 0.24 0.23  -0.71  0.02 1.00      944      779
##  cov[1,1]  1.16   1.05 0.51 0.39   0.60  2.04 1.01     1200      952
##  cov[2,1]  0.91   0.81 0.44 0.33   0.42  1.70 1.01     1011      717
##  cov[1,2]  0.91   0.81 0.44 0.33   0.42  1.70 1.01     1011      717
##  cov[2,2]  1.11   1.01 0.48 0.37   0.58  1.94 1.01     1098      711

distributions

normal dist.

distribution of mean, it's from central limit theorem
y~N(m,s), y(-Inf,Inf), E[y]=m, V[y]=s^2

ex1.stan

data{
  int N;
  vector[N] y;
}
parameters{
  real m;
  real<lower=0> s;
}
model{
  y~normal(m,s);
}
y=rnorm(10,2,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.91    1.16    1.72    2.04    2.40    4.69
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=length(y),y=y)

mdl=cmdstan_model('./ex1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -10.5958 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        7      -6.44154   4.82669e-05   6.98265e-05      0.8803      0.8803       10    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__    -6.44
##      m        2.04
##      s        1.16
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -7.41  -7.07 1.14 0.81 -9.76 -6.32 1.00      755      813
##      m     2.04   2.04 0.48 0.44  1.25  2.83 1.00      752      784
##      s     1.43   1.34 0.42 0.34  0.93  2.20 1.00     1514     1324

Bernoulli dist.

The event with probablity p occur y=1, not occur y=0 
y~Ber(p), y{0,1}, E[y]=p, V[y]=p(1-p)

ex0-0.stan

data{
  int N;
  array[N] int y;
}
parameters{
  real<lower=0,upper=1> p;
}
model{
  p~beta(1,1);
  y~bernoulli(p);
}
data=list(N=9,y=sample(0:1,9,replace=T))

mdl=cmdstan_model('./ex0-0.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -8.44677 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5      -5.72863   0.000161808   6.74689e-07           1           1        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__    -5.73
##      p        0.67
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -7.70  -7.45 0.68 0.32 -8.96 -7.21 1.00     1035     1383
##      p     0.64   0.65 0.13 0.14  0.40  0.84 1.01      880     1045

Poisson dist.

The event occur l times in unit range, the event occur y times. 
y~Po(l), y{0,1...Inf}, E[y]=l, V[y]=l

ex2-2.stan

data{
  int N;
  array[N] int y;
}
parameters{
  real<lower=0> l;
}
model{
  y~poisson(l);
}
generated quantities{
  array[N] int y1;
  for(i in 1:N)
    y1[i]=poisson_rng(l);
}
n=10
y=rpois(n,3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    1.25    2.00    2.90    4.50    7.00
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex2-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -15.7744 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5       1.87661   3.45392e-05   2.78157e-06           1           1        9    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__       1.88
##    l          2.90
##    y1[1]      2.00
##    y1[2]      1.00
##    y1[3]      2.00
##    y1[4]      1.00
##    y1[5]      2.00
##    y1[6]      2.00
##    y1[7]      1.00
##    y1[8]      1.00
##    y1[9]      4.00
##    y1[10]     0.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable mean median   sd  mad   q5  q95 rhat ess_bulk ess_tail
##    lp__   2.46   2.72 0.67 0.33 1.20 2.96 1.00      895     1242
##    l      2.99   2.96 0.54 0.56 2.15 3.92 1.00      470      830
##    y1[1]  2.94   3.00 1.78 1.48 0.00 6.00 1.00     1900     1810
##    y1[2]  2.97   3.00 1.76 1.48 0.00 6.00 1.00     1938     1894
##    y1[3]  2.97   3.00 1.82 1.48 0.00 6.00 1.00     1964     1794
##    y1[4]  3.01   3.00 1.83 1.48 0.00 6.00 1.00     1523     1860
##    y1[5]  2.97   3.00 1.81 1.48 0.00 6.00 1.00     1700     1683
##    y1[6]  2.94   3.00 1.81 1.48 0.00 6.00 1.00     1538     1798
##    y1[7]  3.00   3.00 1.80 1.48 0.00 6.00 1.00     1700     1884
##    y1[8]  2.96   3.00 1.80 1.48 0.00 6.00 1.00     1807     1834
##    y1[9]  2.99   3.00 1.80 1.48 0.00 6.00 1.00     1732     2021
##    y1[10] 2.96   3.00 1.77 1.48 0.00 6.00 1.00     1916     1895

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

binomial dist.

The event with probablity p occur y times on n trials
y~B(n,p), y{0,1...n}, E[y]=np, V[y]=np(1-p)

ex2-3.stan

data{
  int N;
  array[N] int n;
  array[N] int y;
}
parameters{
  real<lower=0> p;
}
model{
  y~binomial(n,p);
}
n=10
n0=floor(runif(n,1,10))
y=rbinom(n,n0,0.3)

data=list(N=n,n=n0,y=y)

mdl=cmdstan_model('./ex2-3.stan')

mle=mdl$optimize(data=data)
## Rejecting initial value: 
##   Error evaluating the log probability at the initial value. 
## Exception: binomial_lpmf: Probability parameter is 3.32476, but must be in the interval [0, 1] (in '/tmp/RtmpXPrRnk/model-8004714516a.stan', line 12, column 2 to column 18) 
## Rejecting initial value: 
##   Error evaluating the log probability at the initial value. 
## Exception: binomial_lpmf: Probability parameter is 1.56971, but must be in the interval [0, 1] (in '/tmp/RtmpXPrRnk/model-8004714516a.stan', line 12, column 2 to column 18) 
## Rejecting initial value: 
##   Error evaluating the log probability at the initial value. 
## Exception: binomial_lpmf: Probability parameter is 2.72465, but must be in the interval [0, 1] (in '/tmp/RtmpXPrRnk/model-8004714516a.stan', line 12, column 2 to column 18) 
## Initial log joint probability = -27.6751 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6      -25.4602   0.000303008   0.000379052      0.9528      0.9528        9    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -25.46
##      p        0.36
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -26.97 -26.71 0.70 0.34 -28.36 -26.46 1.00      965      932
##      p      0.36   0.36 0.07 0.08   0.25   0.49 1.00      673      767

log normal dist.

logalithm of variable follows normal distribution
log y~N(m,s), y>0

ex2-4.stan

data{
  int N;
  vector[N]<lower=0> y;
}
parameters{
  real m;
  real<lower=0> s;
}
model{
  y~lognormal(m,s);
}
n=10
y=rlnorm(n,2,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.9     2.1     4.3     9.8    14.8    34.9
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex2-4.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -37.1816 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       16       -15.534   0.000182535   0.000119359      0.9162      0.9162       20    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -15.53
##      m        1.67
##      s        1.14
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -16.50 -16.17 1.10 0.79 -18.67 -15.44 1.00      804      982
##      m      1.67   1.65 0.46 0.43   0.96   2.47 1.01      465      625
##      s      1.42   1.34 0.41 0.34   0.92   2.17 1.00     1388     1446

exponential dist.

The event occur l times in unit range, the event take time y to occur
y~ex(l), y>0, E[y]=1/l, V[y]=1/l^2

ex2-5.stan

data{
  int N;
  vector[N] y;
}
parameters{
  real<lower=0> l;
}
model{
  y~exponential(l);
}
n=10
y=rexp(n,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.04    0.26    0.92    1.01    1.35    3.19
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex2-5.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -18.0237 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6      -10.0976    0.00104077   9.54535e-05           1           1        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -10.10
##      l        0.99
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -10.55 -10.27 0.68 0.29 -12.00 -10.06 1.00      928     1106
##      l      1.08   1.05 0.32 0.31   0.61   1.65 1.00      687      873

gamma dist.

The event occur l times in unit range, the event take time y to occur a times
y~Ga(a,l), y>0, E[y]=a/l, V[y]=a/l^2

ex2-6.stan

data{
  int N;
  vector[N] y;
}
parameters{
  real<lower=0> a;
  real<lower=0> l;
}
model{
  y~gamma(a,l);
}
n=10
y=rgamma(n,2,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.11    1.54    2.13    2.28    2.65    5.97
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex2-6.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -25.3383 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       10      -17.9853    0.00100604   1.24513e-06           1           1       13    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -17.99
##      a        1.35
##      l        0.59
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -18.98 -18.65 1.06 0.78 -21.17 -17.96 1.01      626      643
##      a      1.66   1.57 0.63 0.61   0.80   2.84 1.01      580      785
##      l      0.78   0.73 0.33 0.30   0.32   1.40 1.01      549      546

beta dist.

use as prior of binomial dist parameter p
p~Be(a,b), p[0,1], E[p]=a/(a+b), V[p]=ab/(a+b)^2/(a+b+1)

ex2-7.stan

data{
  int N;
  vector[N] p;
}
parameters{
  real<lower=0> a;
  real<lower=0> b;
}
model{
  p~beta(a,b);
}
n=20
p=rbeta(n,1,2)
summary(p)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.066   0.163   0.289   0.351   0.441   0.896
par(mfrow=c(1,2))
hist(p,main='p')
density(p) |> plot(main='p')

data=list(N=n,p=p)

mdl=cmdstan_model('./ex2-7.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -132.547 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       11       3.14447   0.000459542   5.50416e-05           1           1       12    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__     3.14
##      a        1.24
##      b        2.13
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable mean median   sd  mad   q5  q95 rhat ess_bulk ess_tail
##      lp__ 3.23   3.54 1.01 0.75 1.18 4.21 1.00     1001     1085
##      a    1.39   1.35 0.38 0.38 0.83 2.07 1.01      577      997
##      b    2.43   2.39 0.69 0.67 1.39 3.64 1.01      604     1025

Dirichlet dist.

use as prior of categorical dist parameter p
p~dir(p0), p[0,1]

ex2-8.stan

data {
  int N;
  int K;
  matrix[N, K] p;
}
parameters {
  simplex[K] p0;
}
model {
  for(i in 1:N){
     p[i]~dirichlet(p0);
  }
}
library(gtools)
n=20
p=rdirichlet(n,c(0.5,0.3,0.2))
summary(p)
##        V1              V2              V3       
##  Min.   :0.002   Min.   :0.000   Min.   :0.000  
##  1st Qu.:0.120   1st Qu.:0.007   1st Qu.:0.006  
##  Median :0.547   Median :0.117   Median :0.033  
##  Mean   :0.542   Mean   :0.328   Mean   :0.130  
##  3rd Qu.:0.930   3rd Qu.:0.673   3rd Qu.:0.166  
##  Max.   :0.998   Max.   :0.998   Max.   :0.896
boxplot(p)

data=list(N=n,K=ncol(p),p=p)

mdl=cmdstan_model('./ex2-8.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = 51.2221 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        4       53.1332    0.00825801    0.00122163           1           1        6    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##     lp__     53.13
##     p0[1]     0.49
##     p0[2]     0.27
##     p0[3]     0.23
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##     lp__  48.70  49.01 0.98 0.71 46.73 49.63 1.00     1054     1197
##     p0[1]  0.48   0.48 0.06 0.06  0.39  0.58 1.00     1965     1405
##     p0[2]  0.28   0.28 0.05 0.05  0.20  0.37 1.00     1412     1167
##     p0[3]  0.24   0.24 0.04 0.04  0.17  0.31 1.00     1726     1408

categorical dist. with parameters following Dirichlet dist.

h~Dir(a), h[0,1], sum(h)=1, a[0,1], sum(a)=1 
y~Cat(h), y[0,1], sum(y)=1

ex2-9.stan

data {
  int N;
  int K;
  array[N] int<lower=1,upper=K> y;
}
parameters {
  simplex[K] a;
  simplex[K] h;
}
model {
  h~dirichlet(a);
  y~categorical(h);
}
library(gtools)
n=20
h=rdirichlet(n,c(0.5,0.3,0.2))
summary(h)
##        V1              V2              V3       
##  Min.   :0.015   Min.   :0.000   Min.   :0.000  
##  1st Qu.:0.126   1st Qu.:0.025   1st Qu.:0.004  
##  Median :0.253   Median :0.214   Median :0.087  
##  Mean   :0.422   Mean   :0.297   Mean   :0.281  
##  3rd Qu.:0.758   3rd Qu.:0.554   3rd Qu.:0.593  
##  Max.   :1.000   Max.   :0.974   Max.   :0.985
c0=c(1,2,3)
for(i in 1:n) y[i]=sample(c0,1,prob=h[i,])
table(y) |> prop.table()
## y
##    1    2    3 
## 0.45 0.35 0.20
data=list(N=n,K=ncol(h),y=y)

mdl=cmdstan_model('./ex2-9.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -28.4628 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       13       -21.589   0.000450354   2.27097e-05           1           1       17    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -21.59
##      a[1]     0.37
##      a[2]     0.34
##      a[3]     0.28
##      h[1]     0.47
##      h[2]     0.35
##      h[3]     0.18
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -30.60 -30.24 1.54 1.38 -33.46 -28.79 1.01      844     1062
##      a[1]   0.34   0.33 0.18 0.19   0.08   0.66 1.00     1379     1261
##      a[2]   0.34   0.33 0.18 0.19   0.08   0.68 1.00     1213      983
##      a[3]   0.31   0.29 0.17 0.18   0.06   0.61 1.00     1061      753
##      h[1]   0.45   0.45 0.11 0.11   0.27   0.63 1.00     3059     1629
##      h[2]   0.35   0.34 0.10 0.11   0.19   0.53 1.00     2371     1592
##      h[3]   0.20   0.19 0.09 0.09   0.08   0.37 1.00     2393     1689